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ABSTRACT 

One  method  of  reducing  truncation  error  in  numerical 
prediction  is  the  use  of  more  accurate  finite-difference 
approximations.   A  relatively  simple  barotropic  model, 
similar  to  the  one  currently  in  use  by  the  Fleet  Numerical 
Weather  Facility  (FNWF),  Monterey,  California,  is  employed 
as  a  basic  program  in  testing  several  higher  order  finite- 
differencing  schemes.   Forecasts  were  computed  to  24  and 
4#  hours  with  modifications  of  the  basic  model,  and 
comparisons  were  made  with  the  analyses  at  verifying  time. 
It  is  the  intent  of  this  study  to  determine  the  effects  of 
some  higher  order  finite-difference  approximations  in  this 
numerical  prediction  model. 

The  author  expresses  his  appreciation  to  Professor 
George  J.  Haltiner  for  his  suggestions,  guidance,  and 
patience.   Appreciation  is  also  expressed  to  the  personnel 
of  FNWF  for  providing  data  and  programming  assistance. 
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I.    Introduction. 

In  any  numerical  method  of  weather  prediction, 
regardless  of  the  simplicity  or  complexity  of  the  forecast 
model,  one  source  of  error  which  cannot  be  entirely 
eliminated  is  the  so-called  truncation  error.   This  error 
is  mathematical  in  nature  and  stems  from  the  use  of  finite- 
difference  approximations  to  the  derivatives  of  functions 
such  as  those  encountered  in  meteorology.   Thompson  [7] 
discusses  the  problem  of  truncation  error  and  outlines 
two  methods  of  reducing  the  magnitude  of  the  error.   The 
first  method  involves  decreasing  the  grid  interval  thereby 
increasing  the  number  of  grid  points  for  a  given  area. 
This  procedure  is  not  considered  in  this  study.   The  second 
method  utilizes  the  "quasi-analytic"  property  of  the 
functions  dealt  with  in  numerical  prediction  which  implies 
that  the  functions  and  their  derivatives  vary  smoothly. 
If  in  reality  these  functions  are  analytic,  uniformly 
convergent  Taylor  series  can  be  written  for  any  such 
functions.   The  inclusion  of  higher  order  terms  of  the 
Taylor  expansion  results  in  more  accurate  finite-difference 
approximations . 

It  is  the  purpose  of  this  study  to  determine  the 
effects  of  some  higher  order  finite-difference  approximations 
in  numerical  weather  prediction  with  the  intention  of 
obtaining  a  more  accurate  prognostication. 


2.   The  Basic  Model 

The  basic  model  used  in  this  study  is  a  simple  model 
described  by  Haltiner  [4]  and  similar  to  the  500-mb 
barotropic-type  forecast  currently  in  operational  use  by 
the  Fleet  Numerical  Weather  Facility  (FNWF),  Monterey, 
California  [2] .   This  model  was  selected  because  of  its 
relative  simplicity,  availability  of  data  in  a  convenient 
form,  and  suitability  in  programming .   It  was  supposed  that 
a  simple  model  would  lend  itself  best  to  this  type  of  study 
In  this  model,  three  empirical  terms  are  incorporated; 
(a)  a  "Helmholtz"  term,  (b)  a  "latitudinal"  term,  and 
(c)  a  "temperature"  term.   Although  these  terms  are 
reported  to  be  relatively  small,  they  are  included  because 
they  have  been  found  to  improve  the  results. 

The  equation  programmed  for  the  basic  model  is  of  the 
form: 

(v*-b)H  =-cJ(z,vh). 

Equation  (1)  is  solved  numerically  by  a  sequential 
relaxation  method.   The  input  parameters  are  the  500-mb 
analysis  and  the  thickness  of  the  layer  between  1000  mb  and 
500  mb,  with  only  the  values  north  of  2N  entering  into  the 
computations.   A  forward  time  step  is  used  initially 
followed  by  centered-diff erences  at  one-hour  intervals.   At 
six-hour  intervals,  slight  smoothing  is  applied  and  the 


ellipticity  criterion  imposed,  after  which  a  forward  time 
step  is  again  initiated  and  the  process  repeated „ 

The  finite-difference  approximations  programmed  in 
the  basic  forecast  scheme  can  be  found  in  most  texts  on 
numerical  methods.   Collatz  [l]  provides  a  comprehensive 
list  of  these  formulas.   Formulas  used  in  the  basic  program 
are  the  5-point  approximation  to  the  Laplacian 


¥"Z0  *  ji   (z,+Z2  +  Z3*Z,-4Z0) 


and  the  5-point  approximation  to  the  Jacobian 


(2) 


J  WO  *  ^R<»i-*.Xvfe>(V*0(v+, 


(3) 


where  the  numerical  subscripts  refer  to  the  grid  position 
with  the  5-point  grid  defined  as 

1.0     2. 


T 

d 

The  time  stepping  formulas  consist  of  the  forward-difference 


z      -  z  +7^1  At 
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and  the   centered-difference 
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t+At 


+   Z  (*1  )   A  t 


't-A-t 
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where  the  subscripts  refer  to  the  time  at  which  the 
quantities  are  computed. 

To  reduce  the  truncation  error  in  the  basic  forecast 
program  many  higher  order  finite-difference  formulas  are 
available.   Only  a  few  of  these  formulas  are  considered 
in  this  investigation.   To  identify  the  modifications  to 
the  basic  program,  letter  designations  were  given  to  the 
changes  which  are  described  in  the  following  section. 

3.   Modification  B. 

The  first  modification  to  the  basic  program  consists 
of  substituting  for  equation  (2)  a  9-point  approximation 
to  the  Laplacian 
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z^^l-^i^^z^^i^y^z^z^. 


(6a) 


Equation  (6a)  can  also  be  written  symbolically  as  suggested 
by  Forsythe  and  Wasow  [3 J  '• 


vz;-^ 


-/ 

Ik    -60  Ho     -/ 
-I 


H. 


(6b) 


The  remainder,  representing  the  residual  between  the 
finite-difference  approximation  and  the  actual  value, 
consists  of  sixth-order  derivatives  and  higher. 


4.   Modification  C. 

The  same  grid  as  presented  in  section  3  can  be  used  to 


derive  a  9-point  approximation  to  the  Jacobian  applying 
the  formula  for  a  5-point  centered-diff erence  along  each 
dimension  of  the  grid.   Modification  C  includes  the 
replacement  of  equation  (3)  by 

jr(*./k)  -  ^[(V^  *+..-**X*i  -»V»iH) 

— 


(7) 


where  the  subscripts  locate  the  positions  on  the  grid  as 
follows : 


S 


a. 


The  remainder  contains  the  product  of  third-order  derivatives 
as  the  leading  term.   In  both  modifications  B  and  C  no  other 
changes  are  made  to  the  basic  program,  and  no  additional 
computer  storage  space  of  any  significance  is  required . 
By  the  inclusion  of  both  9-point  approximations,  equations 
(6)  and  (7),  a  composite  program  is  obtained.   This  dual 
program  is  designated  as  modification  X. 

5.   Modification  D. 

In  the  basic  program,  equations  (4)  and  (5)  are  used 
for  the  first  and  second  hourly  time  steps  respectively. 
For  this  modification,  subsequent  time  steps  are  performed 


by  utilizing  an  approximation  comprised  of  terms  through  the 
third  derivative  as  shown: 

z        _7        ~   tfe)*t+ £.(££)  Utf. 

*-£+A±       £t-At  *-\cbif:/£  3l\£>t*kK      ' 

(d) 

By  the  method  of  divided  differences  described  by  Jennings 

[5]  a  formula  for  equation  (&)  in  terms  of  £  and  ^?can  be 
derived: 

(9) 

Equation  (9)  has  the  same  form  as  the  formula  listed  in 
Milne  [6]  for  a  net  of  four  equally-spaced  points.  The 
remainder  can  be  shown  to  contain  terms  higher  than  the 
third  derivative. 

Note  that  equation  (9)  contains  one  additional  past 
point.  In  programming  this  modification,  one  additional 
field  is  necessary  in  the  computer  memory. 

6.   Modification  E. 

Iterations  after  the  first  two  time  steps  are 
accomplished  in  this  case  by  an  alternative  method  of 
including  the  third-order  term  in  equation  (#).   By 
considering 


£t* 


and  using  a  one-sided  difference  for  the  second  derivative 
in  terms  of  the  first  derivative,  equation  (£)  becomes 

(10) 
Equation  (10)  contains  three  height-change  tendencies  and 
consequently  greater  computer  storage  space  is  required. 
Due  to  limited  computer  storage  space  and  other  considerations 
to  be  pointed  out  in  section  7,  both  modifications  D  and  E 
necessitated  considerable  revision  of  the  program. 

7.    Procedures. 

All  computations  were  performed  with  the  CDC  1604 
all  purpose  computer  using  a  63  x  63  square  grid  with  a  mesh 
length  of  3#1  km.   The  grid  is  superimposed  on  a  polar 
stereographic  projection  of  the  northern  hemisphere  centered 
at  the  North  Pole  and  true  at  60N. 

Imput  data  in  the  form  of  500-mb  and  1000-mb  analyses 
on  magnetic  tape  were  provided  by  FNWF0   The  period  from 
00Z  15  March  1964  through  00Z  20  March  1964  was  used  in 
this  study.   This  period  was  selected  because  the  synoptic 
situation  showed  a  block  positioned  over  northern  Europe, 
a  series  of  waves  moving  rapidly  across  the  Pacific  area, 
and  systems  developing  and  weakening  across  the  North 
American  continent.   This  depicted  a  variety  of  conditions 
for  which  the  programs  could  be  tested. 


The  prognosis  was  carried  out  to  24  and  4#  hours  for 
all  programs  with  the  exception  of  modification  E„ 
Difficulty  was  encountered  with  modification  D  because  of 
overflow  after  36  hours.   By  a  small  amount  of  smoothing 
at  each  time  step  the  program  was  able  to  operate  to  no 
less  than  4S  hours.  Modification  E  overflowed  at  30  hours . 
Smoothing  the  height  field  or  the  tendency  field  was  in- 
effective in  controlling  the  instability,  hence  modification 
E  was  terminated  at  24  hours . 

When  each  series  of  computations  was  completed,  the 
products  were  compared  at  each  grid  point  with  the  analysis 
at  the  appropriate  verifying  times .   The  pillow  and  root 
mean  square  error  were  computed  for  each  prognosis ,      The 
pillow  and  root  mean  square  error  are  denoted  respectively: 

n 


P**;£  (*.-2,)t 


C-l 


and 


where  n  is  the  number  of  grid  points  considered .   For  this 
comparison  all  grid  points  were  used,  hence  n=3969« 

8.        Results  and  Discussion. 

The  results  of  the  comparisons  are  given  in  Table  I. 
It  should  be  noted  that  the  results  are  slightly  biased  in 
favor  of  the  basic  program.   This  is  due  to  the  use  of  the 

d 


TABLE  I 


Pillow  and  RMSE  of  verifying  analysis  and  prognosis 


Prognosis 


Basic 


Mod  B 


Mod  C 


Mod  D 


Mod  E 


Mod  X 


Initial  Time 

24-hr 

P/RMSE 

(meters) 

43-hr 

R/RMSE 

(meters) 

15/00Z 

6.4/31.7 

7.9/50.3 

16/00Z 

2.4/31=1 

3.7/54.0 

17/00Z 

2.1/35.4 

0.3/54.9 

1S/00Z 

-0.3/29.6 

1.3/43.5 

15/00Z 

7.0/47.3 

3. 2/59.7 

16/00Z 

3.1/46.0 

5. 5/77. 3 

15/00Z 

7=6/31.7 

9.4/49.1 

16/00Z 

3.3/30,5 

5.2/52.7 

17/00Z 

3.3/34.1 

4.7/54.0 

lg/OOZ 

0.3/29.9 

3.1/43.1 

15/00Z 

16.8/4S.S 

44.2/97.6 

16/00Z 

12.2/43.3 

39.6/36.0 

15/00Z 

8.3/62.2 

-- 

16/00Z 

4.3/55.2 

— 

15/00Z 

7.3/50.0 

3.3/63.0 

16/00Z 

3.1/47.0 

5.5/74.0 

basic  prognosis  in  the  first  approximation  for  producing 
the  verifying  analysis.   Maps  containing  the  500-mb  height 
field  were  printed  as  a  grid  for  each  set  of  computations 
and  contour  charts  plotted  with  the  CDC  160  and  graph  plotter, 


The  prognoses  from  00Z  16  March  1964  for  each  program, 
typifying  the  products  obtained,  are  given  in  Appendix  I. 
Analyses  of  the  00Z  charts  covering  the  period  concerned  are 
also  presented. 

By  examination  of  the  products  of  the  various  programs 
and  comparison  with  the  analyses  at  verifying  time,  the 
following  observations  were  made: 

(a)  All  programs  show  reasonable  contour  patterns 
with  the  exception  of  modification  D.   It  is  obvious  from 
charts  7  and  £  that  spurious  perturbations  with  a  wave- 
length on  the  order  of  three  grid-intervals  developed  and 
amplified  with  time.   These  waves  appear  to  be  caused  by 
computational  instability  originating  from  the  use  of  the 
finite-difference  formula. 

(b)  The  major  polar  low  was  too  deep  in  all  programs, 
ranging  from  20-60  m  at  24  hours  to  105-160  m  at  4#  hours. 
On  the  other  hand,  the  heights  of  quasi-stationary  systems 
tended  to  be  too  large  by  an  average  of  12  m  for  highs  and 
70  m  for  lows  at  24  hours,  and  43  m  for  highs  and  71  m  for 
lows  at  4#  hours.  The  magnitude  of  the  differences  tended 
to  contribute  to  the  relative  accuracy  of  the  individual 
programs . 

(c)  The  speed  of  moving  systems  was  in  general  slower 
than  the  verifying  analysis  indicated.   Interestingly,  move- 
ments in  modifications  B  and  X  could  hardly  be  detected  in 
most  cases.   Recall  that  modifications  B  and  X  both  contain 
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the  9-point  approximation  to  the  Laplacian.   This  change 
was  suspect  to  possible  programming  error  or  incompatibility 
of  the  modification  and  the  basic  program.   A  check  of  the 
subroutine  for  the  9-point  Laplacian  gave  reasonable 
calculations.   Consequently,  the  latter  reason  is  indicated 
as  the  cause  of  the  results,  especially  in  view  of  the 
tendency  for  spurious  anti-cyclogenesis  to  creep  in.   Anti- 
cyclogenesis  was  controlled  in  the  other  programs. 

(d)   As  expected  for  barotropic  models,  all  programs 
suffered  when  predicting  the  intensification  of  developing 
systems. 

Of  all  the  modifications  tested,  the  only  one  which 
compares  favorably  with  the  basic  program  is  modification 
C  which  utilizes  the  9-point  Jacobian.   There  is  strong 
evidence  that  the  use  of  this  subroutine  can  improve  the 
basic  program.   Due  to  the  limited  computer  time  available, 
the  length  of  each  run  (see  table  II),  and  other  difficulties, 
a  relatively  small  sample  was  obtained.   Because  the  9-point 
Jacobian  program  showed  more  promise,  additional  data  were 
collected  for  modification  C„   Similar  data  for  the  basic 
program  were  also  required  to  determine  if  the  two  methods 
were  comparable.   The  RMSE  for  modification  C  averaged 
slightly  less  than  the  basic  program,  though  the  pillow 
was  somewhat  higher.   There  was  no  significant  difference  in 
running  time. 
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Prognosis 

Basic 

Mod 

B 

Mod 

C 

Mod 

D 

Mod 

E 

Mod 

X 

TABLE  II 
Computer  running  time  for  the  prognoses 

2L-hr   Prog  1+8-hr   Prog 

(minutes)  (minutes ) 

12  24 

14  27 

12  24 

24  47 

27 

16  32 


9.   Conclusions. 

From  this  investigation,  it  can  be  said  that  although 
truncation  error  can  be  reduced  theoretically  through  the 
use  of  higher  order  finite-differences,  the  incorporation 
of  these  approximations  into  a  numerical  prognostication 
scheme  can  lead  to  undesirable  effects  in  some  cases. 
Computational  instability  and  other  difficulties  encountered 
in  this  study  point  out  the  implicit  complexity  of  making 
a  seemingly  simple  change  to  a  forecast  method.   No 
modifications  other  than  the  substitution  of  higher  order 
finite-difference  formulas  were  used,  except  for  the  small 
amount  of  additional  smoothing  previously  mentioned.   By 
adding  empirical  corrections  to  the  model,  it  is  conceivable 
that  the  forecasts  could  be  further  improved,  but  this  was 


12 


not  the  primary  objective  of  this  paper. 

The  limited  results  of  this  study  provide  additional 
avenues  of  research.   One-sided  advection  schemes  and 
numerous  other  stencils  for  the  finite-difference  approxi- 
mations are  suggested  for  further  examination. 
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APPENDIX  I 

The  following  charts  were  produced  during  this  investigation, 
They  are  typical  of  the  prognoses  obtained  by  the  programs 
studied.   Also  included  are  the  analyses  of  the  charts  for 
the  period  covered. 


Charts  1-2 
Charts  3-4 
Charts  5-6 
Charts  7-g 
Chart  9 
Charts  10-11 
Charts  12-17 


24-hr  and  4S-hr  prognoses,  00Z 
16  March  1964  Basic  program 

24-hr  and  4#-hr  prognoses,  00Z 
16  March  1964  Modification  B 

24-hr  and  4#-hr  prognoses,  00Z 
16  March  1964  Modification  C 

24-hr  and  4^-hr  prognoses,  00Z 
16  March  1964  Modification  D 

24-hr  prognoses,  00Z  16  March 
1964  Modification  E 

24-hr  and  4#-hr  prognoses,  00Z 
16  March  1964  Modification  X 

Analyses,  00Z  15  March  1964 
through  00 Z  20  March  1964 
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